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CENTRAL LIMIT THEOREM FOR THE MULTILEVEL 
MONTE CARLO EULER METHOD 

By Mohamed Ben Alaya^ and Ahmed Kebaier^ 

Universite Paris 13 

This paper focuses on studying the multilevel Monte Carlo method 
recently introduced by Giles [Oper. Res. 56 (2008) 607-617] which is 
significantly more efficient than the classical Monte Carlo one. Our 
aim is to prove a central limit theorem of Lindeberg-Feller type for 
the multilevel Monte Carlo method associated with the Euler dis¬ 
cretization scheme. To do so, we prove first a stable law convergence 
theorem, in the spirit of Jacod and Protter [Ann. Probab. 26 (1998) 
267-307], for the Euler scheme error on two consecutive levels of the 
algorithm. This leads to an accurate description of the optimal choice 
of parameters and to an explicit characterization of the limiting vari¬ 
ance in the central limit theorem of the algorithm. A complexity of 
the multilevel Monte Carlo algorithm is carried out. 


1. Introduction. In many applications, in particular in the pricing of fi¬ 
nancial securities, we are interested in the effective computation by Monte 
Carlo methods of the quantity ]E/(A 7 ’), where X := {Xt)o<t<T is a diffusion 
process and / a given function. The Monte Carlo Euler method consists of 
two steps. First, approximate the diffusion process {Xt)o<t<T by the Eu¬ 
ler scheme (A”)o<t<T with time step T/n. Then approximate E/(Ay) by 
where /(A- ) i<i<N is a sample of N independent copies 
of f{Xlf). This approximation is affected, respectively, by a discretization 
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error and a statistical error 

1 ^ 

en-.= E{f{X^)-f{XT)) and - ^J - ]E/(X^). 

i=l 

On one hand, Talay and Tubaro [21] prove that if / is sufficiently smooth, 
then £n ^ cjn with c a given constant and in a more general context, Kebaier 
[17] proves that the rate of convergence of the discretization error can 
be l/n“ for all values of a € [1/2,1] (see, e.g., Kloeden and Platen [18] for 
more details on discretization schemes). On the other hand, the statistical 
error is controlled by the central limit theorem with order 1/^/N. Further, 
the optimal choice of the sample size N in the classical Monte Carlo method 
mainly depends on the order of the discretization error. More precisely, it 
turns out that for = l/n“ the optimal choice of N is n^“. This leads to 
a total complexity in the Monte Carlo method of order Cmc = (see 

Duffie and Glynn [5] for related results). Let us recall that the complexity of 
an algorithm is proportional to the maximum number of basic computations 
performed by this one. Hence, expressing this complexity in terms of the 
discretization error , we get Cmc = En ^ • 

In order to improve the performance of this method, Kebaier introduced a 
two-level Monte Carlo method [17] (called the statistical Romberg method) 
reducing the complexity Cmc while maintaining the convergence of the algo¬ 
rithm. This method uses two Euler schemes with time steps T/n and T/n^, 
/3 € (0,1) and approximates E/(X'r) by 



where X^^ is a second Euler scheme with time step T/rfi and such that the 
Brownian paths used for XJj^ and Xl^f has to be independent of the Brownian 
paths used to simulate . It turns out that for a given discretization error 
£„ = l/n“ (a G [1/2,1]), the optimal choice is obtained for j3 = 1/2, Ni = 
and N 2 = With this choice, the complexity of the statistical 

Romberg method is of order Csr = which is lower 

than the classical complexity in the Monte Carlo method. 

More recently, Giles [8] generalized the statistical Romberg method of 
Kebaier [17] and proposed the multilevel Monte Carlo algorithm, in a similar 
approach to Heinrich’s multilevel method for parametric integration [12] (see 
also Creutzig et ah [3], Dereich [4], Giles [7], Giles, Higham and Mao [9], 
Giles and Szpruch [10], Heinrich [11], Heinrich and Sindambiwe [13] and 
Hutzenthaler, Jentzen and Kloeden [14] for related results). The multilevel 
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Monte Carlo method uses information from a sequence of computations with 
decreasing step sizes and approximates the quantity Kf^Xx) by 


No 


Q- = j^Ef(^h)+E 


k=l 


£=1 


1 


Ni 

k=l 




where the fine discretization step is equal to T/n thereby L = • For 

£ e {1,..., L}, processes {xf^ , xf^ )o<t<T, A: € {1,, Nf], are inde¬ 
pendent copies of (Xj’™ )o<t<T whose components denote the Eu¬ 
ler schemes with time steps m~^T and However, for fixed i, the 

simulation of [xf^ )o<t<T and (X^ )o<t<T has to be based on the same 
Brownian path. Concerning the first empirical mean, processes {X^j^)o<t<T, 
A:e{l,...,Xo}, are independent copies of (X/)o<t<r which denotes the Eu¬ 
ler scheme with time step T. Here, it is important to point out that all 
these E -|- 1 Monte Carlo estimators have to be based on different indepen¬ 
dent samples. Due to the above independence assumption on the paths, the 
variance of the multilevel estimator is given by 

L 

:= Var(Q„) = Var(/(X^)) + 

i=i 

where aj = Var(/(X^™ ) — /(X^™ )). Assuming that the diffusion coeffi¬ 

cients of X and the function / are Lipschitz continuous, then it is easy to 
check, using properties of the Euler scheme that 


(T^ < C2 ^ ^ 

£=0 

for some positive constant C 2 (see Proposition 1 for more details). Giles 
[8] uses this computation in order to find the optimal choice of the multi¬ 
level Monte Carlo parameters. More precisely, to obtain a desired root mean 
squared error (RMSE), say of order l/n“, for the multilevel estimator, Giles 
[8] uses the above computation on cj^ to minimize the total complexity of 
the algorithm. It turns out that the optimal choice is obtained for (see The¬ 
orem 3.1 of [8]) 

(1) Ne = 2 c 2 n^'^(^^ + l]— for ^ e {0,...,E} and E = 

\logm y logm 

Hence, for an error = \ln°^ ^ this optimal choice leads to a complexity 
for the multilevel Monte Garlo Euler method proportional to n^“(logn)^ = 
e“^(logen)^. Interesting numerical tests, comparing three methods (crude 
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Monte Carlo, statistical Romberg and the multilevel Monte Carlo), were 
processed in Korn, Korn and Kroisandt [19]. 

In the present paper, we focus on central limit theorems for the inferred 
error; a question which has not been addressed in previous research. To 
do so, we use techniques adapted to this setting, based on a central limit 
theorem for triangular array (see Theorem 2) together with Toeplitz lemma. 
It is worth to note that our approach improves techniques developed by 
Kebaier [17] in his study of the statistical Romberg method (see Remark 2 
for more details). Hence, our main result is a Lindeberg-Feller central limit 
theorem for the multilevel Monte Carlo Euler algorithm (see Theorem 4). 
Further, this allows us to prove a Berry-Esseen-type bound on our central 
limit theorem. 

In order to show this central limit theorem, we first prove a stable law 
convergence theorem, for the Euler scheme error on two consecutive levels 
rr/~^ and of the type obtained in Jacod and Protter [16]. Indeed, we 
prove the following functional result (see Theorem 3): 


, ^ aS £ ^ OO, 

(m — i)r 

where U is the same limit process given in Theorem 3.2 of Jacod and Protter 
[16]. Our result uses standard tools developed in their paper but it cannot be 
deduced without a specihc and laborious study. Further, their result, namely 



X) jj 


as oo. 


is neither sufficient nor appropriate to prove our Theorem 4, since the mul- 

0 ^ 0 ^ 1 

tilevel Monte Carlo Euler method involves the error process X — X 
rather than X^’^‘ — X. 

Thanks to Theorem 4, we obtain a precise description for the choice of the 
parameters to run the multilevel Monte Carlo Euler method. Afterward, by 
a complexity analysis we obtain the optimal choice for the multilevel Monte 
Carlo Euler method. It turns out that for a total error of order = l/n“ 
the optimal parameters are given by 

(2) Ni= ^"" "^^^ n^^logn for ^ G {0,..., L} and L = 

log m log m 

This leads us to a complexity proportional to n^“(logn)^ = e~‘^{log£n)‘^ 
which is the same order obtained by Giles [8]. By comparing relations (1) 
and (2), we note that our optimal sequence of sample sizes {Ni)o<i<L does 
not depend on any given constant, since our approach is based on proving 
a central limit theorem and not on obtaining an upper bound for the vari¬ 
ance of the algorithm. However, some numerical tests comparing the runtime 
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with respect to the root mean square error, show that we are in line with the 
original work of Giles [8]. Nevertheless, the major advantage of our central 
limit theorem is that it fills the gap in the literature for the multilevel Monte 
Carlo Euler method and allows to construct a more accurate confidence in¬ 
terval compared to the one obtained using Chebyshev’s inequality. All these 
results are stated and proved in Section 3. The next section is devoted to 
recall some useful stochastic limit theorems and to introduce our notation. 

2. General framework. 

2.1. Preliminaries. Let (A„) be a sequence of random variables with 
values in a Polish space E defined on a probability space P). Let 

P) be an extension of (12,7-’,P), and let X be an E-valued random 
variable on the extension. We say that (A„) converges in law to X stably 
and write X^ j£ 

K{Uh{Xn)) ^E{Uh{X)) 

for all h: E —)• M bounded continuous and all bounded random variable U 
on (12,7^). This convergence is obviously stronger than convergence in law 
that we will denote here by According to Section 2 of Jacod [15] and 

Lemma 2.1 of Jacod and Protter [16], we have the following result. 

Lemma 1. Let Vn and V be defined on (12,7^) with values in another 
metric space E'. 

If Vn ^ P, Xn X then {Vn, Xn) (P, X) . 

Conversely, if {V,Xn) {V,X) and V generates the a-field T, we can re¬ 
alize this limit as (P, A) with X defined on an extension of (12,7^, P) and 

^stably 


Now, we recall a result on the convergence of stochastic integrals for¬ 
mulated from Theorem 2.3 in Jacod and Protter [16]. This is a simplified 
version but it is sufficient for our study. Let A” = (A”'’®)i<j<(i be a sequence 
of valued continuous semimartingales with the decomposition 

Xn,i ^ ^ ^ ^ Q < 2 < T, 

where, for each n E N and 1 < * < d, A"’’* is a predictable process with finite 
variation, null at 0 and M"'’* is a martingale null at 0. 

Theorem 1. Assume that the sequence (A”) is such that 

Jo 
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is tight. Let and H be a sequence of adapted, right-continuous and left- 
hand side limited processes all defined on the same filtered probability space. 
If {H,X) then X is a semimartingale with respect to the filtra¬ 

tion generated by the limit process {H,X), and we have (ff”, AC”, J H^dX^) ^ 
{H,XJ HdX). 

We recall also the following Lindeberg-Feller central limit theorem that 
will be used in the sequel (see, e.g., Theorems 7.2 and 7.3 in [1]). 

Theorem 2 (Central limit theorem for triangular array). Let (A:„)nGN 
be a sequence such that —>■ oo as n —)• oo. For each n, let ..., Xn^k„ 
kn independent random variables with finite variance such that E(X„^fc) = 0 
for all k ^ {1,..., kn}- Suppose that the following conditions hold: 

(Al) limn^ooI]fc=iIE|Af„^fcP = o-2,o-> 0. 

(A2) Lindeherg’s condition: for all e > 0, lim„_>.oo x 

l{|W,d>4) = 0- 

^71 

^ AA(0, (T^) asn^oo. 

k=l 

Moreover, if the Xn^k have moments of order p> 2, then the Lindeherg’s 
condition can be obtained by the following one: 

(A3) Lyapunov’s condition: lim„_>.ooE|= 0. 

2.2. The Euler scheme. Let X := {Xt)o<t<T be the process with values 
in solution to 

(3) dXt = h{Xt)dt + a{Xt)dWt, Ao = xGM'=', 

where W = {W^, ..., W^) is a g-dimensional Brownian motion on some given 
filtered probability space B = (At)t>o,P) with (T))t>o is the standard 

Brownian filtration, b and a are, respectively, and valued functions. 
We consider the continuous Euler approximation A” with step 5 = Tjn given 
by 


dXf = biX^„^t))dt + u(A,„(*)) dWi, rjnit) = [t/6]6. 

It is well known that under the global Lipschitz condition 

3Ct > 0, such that, \b{x) — b{y)\ + \cr{x) — cr{y)\ < CT\y — x 
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the Euler scheme satisfies the following property (see, e.g., Bouleau and 
Lepingle [2]): 


Vp> 1, 


{V) 


sup |Xt|, sup \X^\^LP and 

0<t<T 0<t<T 


E 


sup \Xt-X]^\P 
Lo<t<T 


^ Kr>{T) 

— j^p/2 ’ 


Kp{T) > 0. 


Note that according to Theorem 3.1 of Jacod and Protter [16], under the 
weaker condition 


(T~(-b,a) b and a are locally Lipschitz with linear growth, 

we have only the uniform convergence in probability, namely the property 


(V) sup 

0<t<T 

Following the notation of Jacod and Protter [16], we rewrite diffusion (3) as 
follows: 

<? 

dXt = ^{Xt) dYt = Y, ^3 i^t) dY^, 
j=o 

where ipj is the jth column of the matrix a, for 1 < j < q, (po = b and 
Yt := (t, VP/,..., iy/)h Then the continuous Euler approximation X'^ with 
time step 6 = T/n becomes 

(4) dxr = ip{XY))dYt = Y = [t/6]6. 

3=0 


i 

3. The multilevel Monte Carlo Euler method. Let )o<t<T denotes 
the Euler scheme with time step m~^T for G {0,..., L}, where L = logn/ 
logm. Noting that 

L 

(5) Ef{X^)=Ef{X^) + Ynfixf)-f{xf~")), 

i=i 

the multilevel method is to estimate independently by the Monte Carlo 
method each of the expectations on the right-hand side of the above relation. 
Hence, we approximate E/(X^) by 


( 6 ) 


Qn — 


1 

Wo 


No L Nt 

E /(Ma)+E at E(/(4,"/) - fiYj ))■ 

1 _ 1 D _ 1 ^ 1 _ 1 


Here, it is important to point out that all these L + 1 Monte Carlo esti¬ 
mators have to be based on different, independent samples. For each i G 
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the samples ) i<fc<A^ are independent copies of 

) whose components denote the Euler schemes with time 
steps m~^T and and simulated with the same Brownian path. 

Concerning the first empirical mean, the samples fe)i<fc<Ao a-i'e inde¬ 
pendent copies of Xj,. The following result gives us a first description of 
the asymptotic behavior of the variance in the multilevel Monte Carlo Euler 
method. 

Proposition 1. Assume that b and a satisfy condition ('Hb,a)- For a 
Lipschitz continuous function f: —> M, we have 

(7) Var(gj = o('j]X-W^ 

\e=o 


Proof. We have 

L 

Var(QO = Xq-i Var(/(Xi)) + ^Var(/(X^”^') - f{X^f^'~")) 

£=i 

<Xo-^Var(/(X^)) 

L 

+ 2 J]X-HVar(/(X?*^) - /(X^)) + Var(/(X?*^“^) - /(Xt))) 


£=1 


<Xo-^Var(/(X^)) 

L 

+ 2[/]upJ]X-1ie[ sup \xf-Xt\^+ sup \xf-^-Xt\ 


£=1 


'-0<t<T 


0<t<T 


where [/]iip := sup^^.^ ■ We complete the proof by using property 

(P) on the strong convergence of the Euler scheme. □ 


Inequality (7) indicates the dependence of the variance of Qn on the choice 
of the parameters Nq, ... ,Nl. This variance can be smaller than the variance 
of f{Xtf), so that Qn appears as a good candidate for the variance reduction. 

The main result of this section is a Lindeberg-Feller central limit theorem 
(see Theorem 4 below). In order to prove this result, we need to prove first 
a new stable law convergence theorem for the Euler scheme error adapted 
to the setting of multilevel Monte Carlo algorithm. This is crucial and is the 
aim of the following subsection. 
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3.1. Stable convergence. In what follows, we prove a stable law conver¬ 
gence theorem, for the Euler scheme error on two consecutive levels 
and m^, of the type obtained in Jacod and Protter [16]. Our result in Theo¬ 
rem 3 below is an innovative contribution on the Euler scheme error that is 
different and more tricky than the original work by Jacod and Protter [16] 
since it involves the error process ^ rather than — X. 

Note that the study of the error X — X*’™' as £ —>■ oo can be reduced 
to the study of the error X™” — as n oo where X”*” and X” stand 
for the Euler schemes with time steps T/{mn) and T/n constructed on the 
same Brownian path. 

Theorem 3. Assume that b and a are with linear growth then the 
following result holds: 

For allm£N\{0,l} ^st^hiy ^ asn^oo, 

Y (m — IjT 

with {Ut)o<t<T the d-dimensional process satisfying 

f 

where 

(9) — (-^s) Fs,jFs,i with and ‘Ps,i •— 

and {Zt)o<t<T is the valued process solution of the linear equation 

Zt = Id + 'y^ / t€[0,r]. 

1 = 0-^0 

Here, is a dx d matrix with {Xipj)ik is the partial derivative of ipij with 
respect to the kth coordinate, and is a standard -dimensional 

Brownian motion independent ofW. This process is defined on an extension 
(fI,X, (Xt)t>o,R) of the space (0, X, (J'i)t>o,R). 


te[0,T], 


( 8 ) 


Ut = 


V2 


*J=i 


Note that by letting formally m tend to infinity, we recover the Jacod and 
Protter’s result [16]. 


Proof of Theorem 3. Consider the error process , 

defined by 
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Combining relation (4), for both processes X”*” and X”, together with a 
Taylor expansion 

‘‘vr" 

j-o 


where is the dx d matrix whose ith row is the gradient of the real-valued 
function ipij at a point between X” and X™” . Therefore, the equation 

satisfied by C/” can be written as 



In the following, let (Z™"'’”')o<t<r be the valued solution of 



Theorem 48, page 326 in [20], ensures existence of the process ((Zj™'"’’”) ^)o<t<T 
defined as the solution of 



Thanks to Theorem 56, page 333 in the same reference [20], we get 


jjmn,n _ ryvan^n 






dG 


mn^n 

s 


f\z^n.nyl .)"(Y - V” (,,)* 

J 0 _1 


i=i 

<? 


+ 


f {zr-"r' - YTnw) * [■ 

Jo j=l J 


Since the increments of the Euler scheme satisfy 


i=0 
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and 


Y’mn _ ymn _ \ ^ -mn/yri _ yri \ 


i=0 


with ^pl^ = and it is easy to check that 


ur 


( 10 ) 


with 


i,j=l -^0 


n 


t,2 


and 



_ ^ ^mn,n / ^^,,mn,n^y^ _ ^Yi - 

i,j=i ■^o 

9 

ft 

j^mn,n X ^ ymn,n 

t,l / j i 

/ Ki’^-’-{Y:-Yl^^^)ds, 

1=0 

Jo 

9 

ft 

jjmn,n \ ^ rymn,n 

^t,2 — 2^ 

/ i7^~(.-r?„(s))dy/, 

i=i 

^0 

9 

ft _ 

j^mn,n _ ^ ^ ymn,n 

1 Kr^’^iy:-yi^is))ds, 

i=0 J 

0 

9 

ft _ 

Tymn,n \ '• rymn,n 

^t,2 — 2_^ 


9=1 J 

0 


where, for (i,j) G {0,..., g} x {1,..., g}, 
K 


_ (rymn,n\ —i-\ -n -n \ -n \2-n 1 

— \^s ) ^^s,0^s,i ~ 2^(^s,j) ^s,i J) 


ut,T,mn,n / rymn,n\ — ^ -n -n 

tis =[^s ) ‘^s,j^s,u 


and 


mn I 
s.i I ’ 


T,'i,mn,n (rymn,n\ — ^ { -n -mn \ -n \2-' 

V j=i / 

r rymn,n\ — i ,-n ,-zmn 

Now, let us introduce 

Zt = Id+ / 5 ^( 925 ,i dY^)Zs with Lpt,j = Vipj{Xt). 

Jo ,=o 
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Moreover, {{Zt) ^)o<f<r exists and satisfies the following explicit linear 
stochastic differential equation: 

{Zt)-^=h+ j\zs)-^j2{^s,,?ds- fiz^r^j^vs^jdYi. 

J=l i=0 

Thanks to the uniform convergence in probability of the Euler scheme and 
according to Theorem 2.5 in Jacod and Protter [16], we have 

(11) sup |Z™ -Zt\^Q and sup |(Z~)-^ - {Zt)-^\ ^ 0. 

0<t<T 0<t<T 


Furthermore, in relation (10), one can replace, respectively, and 

^i^,j,mn,n common limit given by relation (9). So that relation 

(10) becomes 


( 12 ) 

with 


Ul 


*J=i 


Hld{Y 


'nmn ( 5 ) 


-Kis))dyi+RT 


j^mn,n j^mn^n . T^mn,n . j^mn,n 

+-^^,2 +-^£,3 


jyinn^n 

^£,1 


j^mn,n 

^£,2 


jymn,n 

-^£,3 ’ 


where and i?™"’’”', i e {1,2}, are introduced by relation (10) and 



-i/:-^)(y;-y;^(,))dT/, 


Tjmn,n \ ^ rymn,n 
*0=1 

The remainder term process vanishes with rate y/n in probability. 

More precisely, we have the following convergence result. 


f 




Lemma 2. The rest term introduced in relation (12) is such that 

I /— T-,mn,n\ 

sup \^/nK^ I 
0<t<T 

converges to zero in probability as n tends to infinity. 


For the reader’s convenience, the proof of this lemma is postponed to the 
end of the current subsection. 

The task is now to study the asymptotic behavior of the process given by 
relation (12) 
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In order to study this process, we introduce the martingale process, 

[\yL.is)-ylis))dYi, (z,j) G 
^ 0 

and we proceed to a preliminary calculus of the expectation of its bracket. 
Let {i,j) and G we have: 

• for j //, the bracket = 0 , 

• for j = f and i / i', = 0 , 

• for j = j' and i = i', - ?7n('S)) ds, t G [0,r] and we 

have 

= y^ {r]mn{s) -Vn{s))ds + 0(^-;^j 
™~I |•(mk+£+l)5/m 

= ^ ^ / {Vmnis) - rjnis)) ds + O 

£=0 k=0 d{mk+£)5lm 



(13) 


+ ^ X /I 

■I- k] +0 


E E 

£=0 k=0 

(m — 1)(5^ 
2m 

(m — i)r 
2mn 


m\ m 


[t/6] + O 






t + O 




Having disposed of this preliminary evaluations, we can now study the sta¬ 
ble convergence of ( yj )!<»j<g• virtue of Theorem 2.1 in [15], 

we need to study the asymptotic behavior of both brackets , 

and y/n{M'^’^d for all t G [0,T] and all G {1,..., The 

case j / j' is obvious and we only proceed to prove that: 


• for j = f, ^Y^)t A 0, for all t G [0,r], 

n—)-oo 

• for j = j' and i 7 ^ i', 0, for all t G [0,r], 

n—^oo 

• for j = j' and i = i', ^ ^ [ 0 ,r]. 

n—>-oo 

For the first point, we consider the convergence 
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= 2 


/ g{s,u) 

J0<s<u<t 


ds du 


with 

dmni'^') 'dmni^s') A Tjniu) 

(14) 

- T7„(s) A rymn(tt) + ??n(s) A ??„('«)• 

It is worthy to note that 

(15) gnis) < r]mnis) < s < IJniu) < r]mniu) < U Vs<r?„(ri). 

Hence, g{s,u) = 0, for s < T]n{u), g{s,u) = gmnis) - Vnis), for gn{u) <s<u, 
and 


E(M’"’fo,y^)^^ = 2 


' 0<rin{u)<s<u<t 


iVmnis) - gn{s))dsdu 


T /■* 

<2— / (u — gniu)) du 

n Jo 


<2^t. 


This yields the desired result. Concerning the second point, the norm is 
given by 




l(ylAs)-Kis))KAs)-<is))ds^ 
/ g{s,u)'^ ds du, 

J0<s<u<t 


= 2 

10<s<u<t 

with the same function g given in relation (14). Using the properties of 
function g developed above, we have in the same manner 


E(M”’fo,M"’* ’Tt=2 


I 

Jo< 


t-^-i {gmn{s)-Vnis)) dsdu<2—t, 

0<r]„{u)<s<u<t ^ 


which proves our claim. For the last point, that is the essential one, we use 
the development of given by relation (13) to get 

2 


E( n(M^’fo)^- 


(16) 


2m 


= „^E(M”-M)f-bi^yT!(= + o(i). 
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Otherwise, we have 




(17) 


with 

(18) 


/ / iE((y, 

Jo Jo 

Jo< 


I _ y* 'i 




-Y1 


Vmn (u) rjn («) 


I ds du 


= 2 / h{s, u) ds du 

/ 0 <S<-!i<t 


h(,,u) = iE((yy,,, - - V„(.)) )■ 

On one hand, for s < r]n{u), using property (15) together with the indepen¬ 
dence of the increments r \ — Y^ , . and Y^ , .—Y^ , yields 

rimn(s) Vn(s) rjmn(U) r)n.(u)’ J 

h{s,u) = {r]mn{s) - r]nis)){r]mniu) -T]niu)). 

On the other hand, in relation (18) we use the Cauchy-Schwarz inequality 
to get h{s,u) = 0{^) and this yields 

1 


/ h{s,u) ds du = 0\ 

j 0<rjn(u)<s<u<t 

Now, noting that {r]mnis) - 'nn{s)){r]mniu) - r]niu)) = 0{:^), relation (17) 
becomes 

E((M"'’*’-^')j) = 2 f (rimnis)-rin{s)){'ijmn{u)-r]n{u))dsdu + o(-^ 

Jo<s<u<t 


id mn (•s) -rinis))ds 



Once again thanks to the development of E((M”’*’'^) 4 ) given by relation (13), 
we deduce that 


(19) 




(m- l)2r2 


t‘^ + 0 



By (16) and (19), we deduce the convergence in of n{M^’'^d)t toward 
Theorem 2.1 in Jacod [15], converges 

in law stably to a standard g^-dimensional Brownian motion 
independent of W. Consequently, by Lemma 1 and Theorem 1, we obtain 










16 


M. BEN ALAYA AND A. KEBAIER 


Finally, we complete the proof using relations (11), (12), Lemma 2 and once 
again Lemma 1 to obtain 


mn 


[m — 1)T 


jjmn,n > stably jj 


where = dBl^. 

V2 Jo □ 


Proof of Lemma 2. At first, we prove the uniform convergence in 
probability toward zero of the normalized rest terms for z G {1,2}. 

The convergence of y/nP™”’” i G {1,2} is a straightforward consequence of 
the previous one. The main part of these rest terms can be represented 
as integrals with respect to three types of supermartingales that can be 
classified through the following three cases: 

= {s-i]n{s))ds, 

Jo 

= Vn [ is-T]nis)) dYj, 

Jo 

where (z,j) G{l,...,g}^ and t G [0,T]. In the first case, the supermartingale 
is deterministic of finite variation and its total variation on the interval [0, T] 
has the following expression: 

[ \dD^’^’^\ = ^/n [ {s-r]nis))ds<^. 

Jo Jo Vn 

So, the process converges to 0 and is tight. In the second case, for 

i G {l,...,g}, the supermartingale is also of finite variation and its total 
variation on the interval [0,T] has the following expression: 

£\dDf -°\*. 

It is clear that sup„E(/p^ | dDs’^’^\) < oo, which ensures the tightness of the 
process Therefore, we only need to establish the convergence of 

toward 0 in L^(ll), for t G [0,T]. In fact, we have 

E{{Dr’^f) = 2n [ E{(Y:-Yl^^^){Y:-Yl^^^))dsdu. 

J0<S<U<t 

When s < r]n{u), we have ry„(s) < s < T]n{u) < u and by independence of the 
Brownian motion increments, we deduce that the integrand term is equal to 




THE MULTILEVEL MONTE CARLO METHOD 


17 


0. Otherwise, when s > T]n{u), we apply the Cauchy-Schwarz inequality to 
get 

g((^n,i,0)2) ^ 22" f (^u — rin{u)) du < 2 — t. 

Jo n 

It follows from all these that ^ q. In the last case, for j G {1,..., q}, 

the process M”’ is a square integrable martingale and its bracket has the 
following expression: 

= n /" is-r]nis)fds<—. 

Jo n 

It is clear that sup„E(M'^’^’-^)'r < oo, so we deduce the tightness of the 
process ^^d the convergence ^ g. 

Now thanks to property (V) and relation (11), it is easy to check that the 
integrand processes and }j^d,mn,n^ introduced in relation (10), con¬ 

verge uniformly in probability to their respective limits Kl = {Zs)~^{(ps,o^s,i — 
and = {Zs)~^ips,j(ps,i-, where ipsj = ^^j{Xs) and (ps,i = 
(pi{Xs). Therefore, by Theorem 1 we deduce that the integral processes given 
by 



Kl 


\yJ-YJi^{s))ds and 


jjOJ,mn,n 


{s - rin{s)) dYj 


vanish. Consequently, we conclude using relation (11) that y/nR^^’^ ^ 0 for 

*€{ 1 , 2 }. 

We now proceed to prove that ^ 0. The convergence of the process 

j^mn,n Q jg obviously obtained from the previous one. The main part 

of this rest term can be represented as a stochastic integral with respect to 
the martingale process given by 


with {i,j) G {1,..., gj X {1,..., gj. It was proven in Jacod and Protter [16] 
that 


^stably 


Tf’ 


where (!?*■’)i<ij<q is a standard g^-dimensional Brownian motion defined on 
an extension probability space (0,^'^ of (fljT', (T})t>o,P), which 

is independent of W. Thanks to property (V) and relation (11), the integrand 
process _ 2 /'*J g and once again by Theorem 1 we deduce that 
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the integral processes given by 



vanish. All this allows us to conclude using relation (11). □ 


3.2. Central limit theorem. Let us recall that the multilevel Monte Carlo 
method uses information from a sequence of computations with decreasing 
step sizes and approximates the quantity E/(A 7 ’) by 



log) Tl 

m e N \ {0,1} and L = --. 

logm 

In the same way as in the case of a crude Monte Carlo estimation, let us 
assume that the discretization error 


en = Ef{X^)-Ef{XT) 


is of order 1/n" for any a G [1/2,1]. Taking advantage from the limit the¬ 
orem proven in the above section, we are now able to establish a central 
limit theorem of Lindeberg-Feller type on the multilevel Monte Carlo Euler 
method. To do so, we introduce a real sequence of positive terms 

such that 


(W) 


L 

lim > Of = oo and 

L^oo ^ 

£=1 


lim 




=0 


for p > 2 


and we assume that the sample size depends on the rest of parameters 
by the relation 


( 20 ) 


Ni = 


n^^{m — 1)T 




■E 

i=i 


a£, 


^ G {0,..., L} and L 


log re 
logrre 


We choose this form for Ni because it is a generic form allowing us a straight¬ 
forward use of Toeplitz lemma that is a crucial tool used in the proof of our 
central limit theorem. Indeed, property (W) implies that if is a se¬ 

quence converging to x G M as tends to infinity then 


lim 

L^+oo Y!,i=i ai 


= X. 







THE MULTILEVEL MONTE CARLO METHOD 


19 


In the sequel, we will denote by E, respectively, Var the expectation, respec¬ 
tively, the variance defined on the probability space (O, P) introduced in 
Theorem 3. We can now state the central limit theorem under strengthened 
conditions on the diffusion coefficients. 

Theorem 4. Assume that b and a are functions satisfying the global 
Lipschitz condition Let f be a real-valued function satisfying 

i'Hf) \f{x)-f{y)\<C{l + \x\P + \y\P)\x-y\ for some C,p > 0. 

Assume '¥‘{Xt ^ Vj) = 0, where Vj := {x G f is differentiable at x}, and 
that for some a G [1/2,1] we have 

{T-Le ^) lim n°'£n = Cf (T, a). 

n—>-cx) 

Then, for the choice of Ni, £ G {0,1,..., L} given by equation (20), we have 

- E(/(Xt))) ^ W(C/(r,a),cT2) 

with cj^ = Var(V/(X 7 ’).t/'r) and M{Cf{T,a),a‘^) denotes a normal distribu¬ 
tion. 


The global Lipschitz condition {TLb,a) seems to be essential to establish our 
result, since it ensures property {V). Otherwise, Hutzenthaler, Jentzen and 
Kloeden [14] prove that under weaker conditions on b and a the multilevel 
Monte Carlo Euler method may diverges whereas the crude Monte Carlo 
method converges. 


Proof of Theorem 4. To simplify our notation, we give the proof 
for a = 1, the case a G [1/2,1) is a straightforward deduction. Combining 
relations (5) and (6) together, we get 

Q^-E{f{XT)) = Qi + Ql + en, 

where 

. No 

k=l 

L Ni 

Qn = Y.N- - /(X^™'"))). 

e=i ^ ^ k=i 

Using assumption {Tie^), we obviously obtain the term Cf(T, a) in the limit. 
Taking A^o = ^ apply the classical central limit the¬ 

orem to Qli- Then we have nQ\ -A 0. Finally, we have only to study the 
convergence of nQ^ and we will conclude by establishing 
nQl^Af{0,^T{Vf{XT).UT)). 
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To do so, we plan to use Theorem 2 with the Lyapunov condition and we 
set 


V _^ \ ' ryvnf-,r 


Ne 


and 


( 21 ) 


k=l 


z^‘r‘-‘ := /(4f) - - e(/( 4'™') - /(Xi;’"'")). 


^T,k ■— J y^'-T,k ! J y^^T,k ! 

In other words, we will check the following conditions: 


lim„^oo Eii nXn,if = Var(V/(XT).C/T). 

(Lyapunov condition) there exists p> 2 such that lim„_).oo Yli^=i = 

0 . 

For the first one, we have 
L L 

^E(X„,^)2 = ^Var(X„,,) 


^=l 


( 22 ) 


£=i 

^ 2 
n 


EiVar(Z- 


^=l 


1 ^ e. 

1 ^ 


Var(Z” ). 


Otherwise, since P(X'r ^ 'Df) = 0, applying the Taylor expansion theorem 
twice we get 






= V/(Xt).17™'’'”^ ‘ + - XT)e{XT,X^rf^^ - Xt) 

- -XT)e{XT,xff^'~' -Xt). 

The function e is given by the Taylor-Young expansion, so it satisfies 
e(Y'r,x£’™' — Xr)-'^0 and e{XT,X^^ — Xt) -^0. By property {V), 

£—>■00 £—^00 

we get the tightness of ^ (m-i)T - Xt) and ^ - Xt) 

and then we deduce 


£,m^ 




[m — 1)T 


((Y^’”^^ - Xr)e(YT, - Xt) 


£,mr 


- (Y^’™^ ‘ - Yr)e(Yr, Y^’™' ‘ - Y^)) ^0. 




£—>■00 
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So, according to Lemma 1 and Theorem 3 we conclude that 




(23) 


[m — 1)T 


as oo. 


Using (T-Lf) it follows from property (V) that 
Ve>0 


supE 

£ 


m’- 


(m — 1)T 
We deduce using relation (23) that 


(/(X^’”^V/( 4 ™'"')) 


2 + 6 ’ 


< oo. 


E 


(m — i)r 
Consequently, 


(/(4’"^V/(4’”^'''))V^IE(V/(^t).Ut)"<oo 


for /c G {1,2}. 


m 


- -— Var(Z” ) ^ Var(V/(XT).C/T) < oo. 

[m — 1)1 

Hence, combining this result together with relation (22), we obtain the first 
condition using Toeplitz lemma. Concerning the second one, by Burkholder’s 
inequality and elementary computations, we get for p> 2 


T)^ 

(24) E|X„,,|P = ^E 


Ni 

^T,l 

£=l 


< n _ _ 

- 




where Cp is a numerical constant depending only on p. Otherwise, property 
{V) ensures the existence of a constant Kp > 0 such that 


E\z: 


^ IP 


T,1 


Therefore, 


(25) 


^E\Xr,,i\P<Cp^ 


m 




p £/2 ■ 


e=i 






Cr, 


(Z)r=i i=i 


Taf 0. 

^ n—^rxD 


This completes the proof. □ 
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Remark 1. From Theorem 2, page 544 in [6], we prove a Berry-Esseen- 
type bound on our central limit theorem. This improves the relevance of the 
above result. Indeed, take a = 1 as in the proof, for Xn^ = nQj^ and Xn/ 
given by relation (21), with £ € {1,..., T}, put 


4 = 




Pn — 




n,r| 


£=0 


^=0 


and denote by the distribution function of n{Qn — E/(X^))/s„. Then for 
all X G M and n G N* 


(26) 


\Fn{x) - G{x)\ < 6 


Pn 


where G is the distribution function of a standard Gaussian random variable. 
If we interpret the output of the above inequality as sum of independent 
individual path simulation, we get 

1 




(”i- i)rEr=i«^ 


X l^ao Var(/(Xi)) + E W Var(/(X^”^ ) - f{xff^ ~ )) j . 

According to the above proof, it is clear that Sn behaves like a constant 
but getting lower bounds for Sn seems not to be a common result to our 
knowledge. Concerning taking p = 3 in both inequalities (24) and (25) 
gives us an upper bound. In fact, when / is Lipschitz, there exists a positive 
constant G depending on 6, a, T and / such that 




3/2 


For the optimal choice ag = l, given in the below subsection, the obtained 
Berry-Esseen-type bound is of order 1 / yjlogn. 


Remark 2. Note that the above proof differs from the ones in Kebaier 
[17]. In fact, here our proof is based on the central limit theorem for triangu¬ 
lar array which is adapted to the form of the multilevel estimator, whereas 
Kebaier used another approach based on studying the associated charac¬ 
teristic function. Further, this latter approach needs a control on the third 
moment, whereas we only need to control a moment strictly greater than 
two. Also, it is worth to note that the limit variance in Theorem 4 is smaller 
than the limit variance in Theorem 3.2 obtained by Kebaier in [17]. 
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3.3. Complexity analysis. From a complexity analysis point of view, we 
can interpret Theorem 4 as follows. For a total error of order 1/n", the com¬ 
putational effort necessary to run the multilevel Monte Carlo Euler method 
is given by the sequence of sample sizes specified by relation (20). The as¬ 
sociated time complexity is given by 

Cmmc = C X -I- ^ N£{m^ + with C > 0 


= Cx 


/ — 1)T 

Y ao 


L 

'^ai + 


£=1 


m 



The minimum of the second term of this complexity is reached for the choice 
of weights a| = 1, ^ G {1,..., L}, since the Cauchy-Schwarz inequality en¬ 
sures that ^ O'l: and the optimal complexity for the multi¬ 

level Monte Carlo Euler method is given by 


Cmmc = C X 


{m-l)T 2 q,, , (m^ - l)r 2 a/, x 2 

n logn -I- ——n ilogn) 


oo log m 
= 0(n^“(logn)^). 


m(logm)2 


It turns out that for a given discretization error £^ = 1 / 71 “ to be achieved 
the complexity is given by Cmmc = O{s~'^{log Sn)"^)- Note that this optimal 
choice = 1, ^ G {1,..., L}, with taking oq = 1 corresponds to the sample 
sizes given by 

iV,= ("'~l)^ n2Mogn, £g{0,...,L}. 

log m 

Hence, our optimal choice is consistent with that proposed by Giles [8]. Nev¬ 
ertheless, unlike the parameters obtained by Giles [8] for the same setting 
[see relation (1)], our optimal choice of the sample sizes N£,i G {1,..., L} 
does not depend on any given constant, since our approach is based on prov¬ 
ing a central limit theorem and not on getting upper bounds for the variance. 
Otherwise, for the same error of order = 1/n" the optimal complexity of 
a Monte Carlo method is given by 

CMC=0(n2“+i) = 0(e-2-V«) 


which is clearly larger than Cmmc- So, we deduce that the multilevel method 
is more efficient. Also, note that the optimal choice of the parameter m is 
obtained for m* = 7. Otherwise, any choice Nq = n^"(logn)^, 0 < (3 < 2, 
leads to the same result. Some numerical tests comparing original Giles 
work [8] with the one of us show that both error rates are in line. Here in 
Eigure 1, we make a simple log-log scale plot of CPU time with respect to 
the root mean square error, for European call and with Nq = n^“(logn)^'®. 
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RMSE 


Fig. 1. Comparison of both routines. 


It is worth to note that the advantage of the central limit theorem is to 
construct a more accurate confidence interval. In fact, for a given root mean 
square error RMSE, the radius of the 90%-confidence interval constructed 
by the central limit theorem is 1.64 x RMSE. However, without this latter 
result, one can only use Chebyshev’s inequality which yields a radius equal 
to 3.16 X RMSE. Finally, note that, taking a = 1/2 still gives the optimal 
rate and allows us to cancel the bias in the central limit theorem due to the 
Euler discretization. 

4. Conclusion. The multilevel Monte Carlo algorithm is a method that 
can be used in a general framework: as soon as we use a discretization scheme 
in order to compute quantities such as we can implement the statis¬ 

tical multilevel algorithm. And this is worth because it is an efficient method 
according to the original work by Giles [8]. The central limit theorems de¬ 
rived in this paper fill the gap in literature and confirm superiority of the 
multilevel method over the classical Monte Carlo approach. 
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